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1.0  INTRODUCTION 

Chip-scale  multifunctional  photonic  integrated  circuits  (PICs)  arc  critical  for  realizing 
high  performance  components  required  in  modem  optical  data  transmission  links  and 
telecommunication  systems.  Such  PICs  are  expected  to  play  a  crucial  role  in  next  generation 
smart  and  adaptive  sensors  that  would  require  the  development  of  high  performance  optical 
signal  processors  and  data  routing  technologies  that  will  facilitate  efficient  and  fast  transfer  of 
information,  on-demand  connection  and  ultrafast  data  processing.  They  are  also  essential  for 
reducing  the  production  cost  of  these  systems.  An  important  sub-component  in  realizing  practical 
all-optieal  signal  processor  is  an  optical  memory  (buffer)  element.  Most  suggestions  for  all- 
optieal  memories  and  flip-flops  are  based  on  the  phenomena  of  optical  bistability,  when  a 
nonlinear  optical  system  can  operate  at  two  different  output  intensities  for  the  same  input 
intensity.  In  most  cases  the  nonlinearity  in  semiconductor  based  structures  is  due  to  an  intensity 
dependent  change  of  refractive  index  produced  by  photo-generated  electron-hall  pairs.  Switching 
times  of  these  devices  is  limited  by  recombination  times  of  the  carriers,  which  is  of  the  order  of 
nanoseconds. 

In  the  Phase  1  efforts  undertaken  at  Hybrid  Photonics  in  cooperation  with  Queens  College 
of  CUNY  we  investigated  optical  bistable  devices  based  on  novel  physical  principles,  which 
distinguish  our  efforts  from  all  previously  suggested  approaches.  Two  configurations  were 
studied:  (i)  two  active  waveguides  side  coupled  to  a  ring  resonator,  and  (ii)  two  coupled  active 
microdisk  resonators.  In  both  configurations  CdSe  core-shell  quantum  dots  (QDs)  were  used  as 
active  gain  media. 

In  Section  2,  we  discuss  our  theoretical  efforts  at  understanding  the  underlying  physical 
mechanism  for  the  bistability  and  the  modeling  efforts  to  design  the  bistable  device.  Following 
this  in  Section  3,  we  describe  the  experimental  work  on  the  two  device  configurations  and  the 
results.  Finally,  in  Section  4,  we  summarize  our  findings  and  discuss  the  future  plans. 

2.0  THEORETICAL  MODEL  OF  BISTABILITY  IN  THE  COUPLED  WAVEGUIDE¬ 
RING  GEOMETRY 

The  first  step  in  designing  an  efficient  device  demonstrating  the  bistable  behavior 
observed  in  preliminary  experiments  is  to  achieve  understanding  of  underlying  physical  causes 
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of  the  observed  behavior  and  to  develop  a  computer  model  allowing  us  to  relate  teehnieal 
eharaeteristies  of  the  deviee  to  its  structural  parameters.  This  step  is  divided  in  two  subtasks:  1) 
to  develop  a  theoretical  model  of  the  modes  of  the  waveguides-resonator  system,  and  2)  to 
develop  a  dynamic  model  of  lasing  in  this  structure. 

2.1.  Modes  and  resonant  frequencies  of  the  system 


The  role  of  the  laser  “cavity’’  in  our  system  is  played  by  two  waveguides  coupled  to  each 
other  through  the  ring  resonator.  The  feedback  is  provided  by  reflection  of  light  from  the 
eleaved  edges  of  the  waveguides  and  redirection  of  the  reflected  light  to  the  second  waveguide 
by  coupling  to  the  ring  resonator.  We  present  the  field  in  each  waveguide  as 

Ea{z)  =  Aa[rae<!ye-^],z<La  (1) 

where  Aa  and  qa  are  stationary  amplitudes  and  propagation  constants  of  the  field  in  the 
waveguides  (a  =  1,2),  2  is  the  coordinate  along  each  of  the  waveguide  with  zero  chosen  at  the 
eleaved  edges,  and  La  are  distances  from  the  edges  of  the  respective  waveguides  to  the  coupling 


region  between  the  waveguides  and  the  ring  (sinec  the  size  of  the  eoupler  is  mueh  smaller  than 
the  wavelength,  we  ean  assume  that  the  coupling  oeeurs  at  a  single  point),  and  ra  are  reflection 

coefficients  from  the  eleaved  edges  of  the  waveguides.  Using  standard  model  of  the  coupler  we 
derived  the  relation  between  the  amplitudes  of  the  field  in  the  1st  and  2nd  waveguides  in  the  form: 

A,  =rj(n)rlA/°,''4:kl 


where 


*(«)- 


p(Q)  = 


\K\2  A 

l-|/|2(l-^yo,n>-^2 

\k\2 

l-|/|2(l-£2)e-i*,n)-^'2 


(3) 


In  these  expressions  we  introduced  the  phase  ehange  of  the  field  inside  the  ring 


FlnL  Fl  lire 

O  =  m,  — — -  =  lirm,  — ;  Qr  = - 

2  c  2  FI  n  L 


(4) 
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where  nr  and  Lr  are  refractive  index  and  the  circumference  of  the  ring  resonator,  c  is  the  speed 
of  light,  and  m2  is  an  arbitrary  integer.  We  also  introduced  ring’s  absorption  coefficient,  yr , 

parameters  of  the  coupler  t,K,  which  obey  energy  conservation  condition |/  '+|/q  =  1 ,  and 
losses  at  the  coupler  .  Combining  Equation  (2)  and  (3)  we  obtain  for  real  and  imaginary  parts 
of  the  propagation  constants  qx  ,  =  q\  2  +  iq[  2 : 


t/'L,  +q'2L2  =  nmx 
q'A+q'2L2=  UnS(Q) 


(5) 


where 


y  L 

/  r  r 


S  = 


fc-f)' 


\  +  R2\t\4e 


-aL 


-2(l  -^2)|/p  e  ‘  r^r  cos 


2k  — 


-  n 


(6) 


r  J 


Equation  (5)  shows  that  the  resonance  frequencies  of  the  effective  “cavity”  in  this  configuration 
are  determined  by  properties  of  the  waveguides,  while  the  ring  resonator  affects  only  the 
effective  losses  of  the  system.  The  normalized  modes  of  the  cavity  arc  defined  as 


Ua(z)  =  - 


+e"^: 


(7) 


Ik- 


*•*  +e-'q° 


2.2.  Dynamic  equations 

With  normal  modes  and  resonant  frequencies  determined  we  can  derive  dynamic 
equations  for  field  amplitudes  in  each  of  the  waveguides.  Wc  present  field  in  each  waveguide  as 

Ea{t)  =  \[A^t)e  i0‘^(,)Ua(z)+c.c]  (8) 

where  a  =  1,2,  /loand  y  are  slow  changing  amplitude  and  the  phase  of  the  field,  Qis  an 

unknown  lasing  frequency,  and  use  standard  semielassieal  lasing  theory  [1]  to  derive  dynamic 
equations  for  the  amplitude  and  the  phase  in  each  waveguide: 
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(9) 


where  v'aand  v'are  real  and  imaginary  parts  of  the  resonance  frequencies  of  the  lasing  cavity,  g 
is  normalized  oscillator  strength  of  optical  transitions  responsible  for  lasing,  g(Q)is  the  gain 
profile  of  the  active  medium  centered  at  the  transition  atomic  frequency  ry0with  effective  width 
yL  ,  A p0  is  unsaturated  population  inversion,  which  characterizes  the  pumping  intensity,  Rs  is 
saturation  parameter,  and 


(10) 


is  responsible  for  the  gain  saturation.  For  the  parameters  of  the  structure  realized  experimentally 


the  lasing  is  expected  to  occur  at  a  frequency  far  above  the  threshold  frequency  of  the  waveguide 


so  that  we  can  assume  linear  relationship  between  the  propagation  constants  of  the  waveguide 
and  the  respective  frequencies:  va  ss  cKqa  ,  where  c„  is  the  speed  of  light  in  the  waveguide.  In  this 

case  we  can  combine  Equation  (9)  with  Equation  (5)  to  obtain  closed  system  of  dynamic 
equations  in  the  form 


=  -- ^lnS(fi) 
2  L  y  ' 


where  we  assumed  that  the  waveguides  are  identical.  In  the  stationary  regime  these  equations 
give  an  equation  for  the  lasing  frequency: 


(12) 


Owning  to  the  oscillating  nature  of  the  ring  resonance  parameter  5  this  equation  might 
admits  multiple  solutions  (see  Figure  1  below).  We  have  developed  a  computer  code  written  in 
Matlab  for  solving  Equation  (12)  and  determined  that  for  typical  parameters  of  the  structure  this 
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equation  indeed  might  have  three  solutions  within  the  gain  width  of  the  atomic  transition.  Since 
the  gam  parameter  for  different  frequencies  is  different  each  of  the  possible  frequencies  will 
result  in  different  lasing  intensities.  Thus,  the  presence  of  multiple  solutions  for  the  lasing 
frequency  is  responsible  for  the  multistable  behavior  of  lasing  output  in  the  structure  under 
consideration.  In  order  to  determine  stability  of  possible  solutions  we  need  to  analyze  the 
dynamic  behavior  of  the  amplitudes  in  the  vicinity  of  possible  stationary  values  of  the 
amplitudes.  To  this  end  we  use  the  relation  between  field  amplitudes  in  different  waveguides 
given  by  Equation  (2)  in  order  to  derive  a  dynamic  equation  for  a  single  intensity,  which  can  be 
presented  in  the  following  form 


I  +  /3‘lg‘\U(zf  g(£2)  1  I  /?I|(/(r)|!*(n) 


dz+  — lnS’(Q) 
n 


(13) 


where  /  is  dimensionless  intensity,  and  (3  =  e 


Frequency  is  also  a  dynamic  variable  whose 


initial  dynamics  occurs  over  very  short  times  of  the  order  of  y±'  and  can  be  neglected.  In  this 
case  the  frequency  in  Equation(13)  can  be  related  to  the  instantaneous  value  of  the  intensity  as 

1  \U(z)f  \Ujzf 


mcoB  =  Cl  < 


aaJ| 


1  +  \U(z)\  g(Q)  1  +  Ig2  |t/(z)f  g(Q) 


dz 


(14) 


and  may  evolve  to  one  of  the  possible  stationary  values  depending  on  their  stability  and  the 
initial  conditions.  To  study  the  stability  of  possible  stationary  solutions  we  used  Equations  (13) 
and  (14)  in  the  small  signal  and  plot  the  right  hand  side  of  Equation  (13)  as  a  function  of 
intensity  (Fig.  1).  It  passes  through  zero  at  the  stationary  solutions  and  the  sign  of  the  slopes  of 
this  function  in  the  vicinities  of  zeroes  indicates  the  stability  of  the  respective  solution:  negative 
slope  corresponds  to  stable  solutions,  and  positive  to  the  unstable  ones.  The  black  curve  in  this 
Figure  corresponds  to  a  situation  when  a  zero  intensity  (nonlasing)  solution  is  unstable,  and  there 
is  a  single  stable  lasing  solution  with  nonzero  intensity.  With  increase  of  pumping  (blue  curve) 
the  stability  of  the  solutions  changes:  zero  intensity  (non-lasing)  solution  becomes  stable  again 
and  coexists  with  a  stable  lasing  solution.  This  is  a  rather  unusual  situation,  when  depending  on 
initial  conditions  the  system  can  either  laze  or  not  at  the  same  pumping  level.  Further  increase  of 
pumping  changes  the  situation  again.  Now  the  non-lasing  solution  becomes  unstable  and  two 
stable  lasing  solutions  emerge.  In  order  to  illustrate  this  situation  more  clearly  we  plot  in  Figure 
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2  the  time  dependence  of  intensity  for  the  same  three  pumping  intensities  as  in  Figure  1,  and  for 
small  initial  intensity,  /0  =0.1. 


Figure  1:  Time  Derivative  of  Intensity  as  Function  of  Intensity.  Zeroes  Show  Stationary  Solutions 


f  igure  2:  Time  Dependence  of  Intensity  for  Small  Initial  Intensity 


One  can  see  that  for  small  pumping  (black  curve)  the  intensity  approaches  a  non-zero 
value,  for  the  second  level  of  pumping  (blue  curve)  the  lasing  is  quenched  for  this  initial 
condition,  while  the  further  increase  of  pumping  results  again  in  a  non-zero  stationary  intensity. 
In  order  to  demonstrate  the  existence  of  the  second  stable  solution  at  blue  pumping,  we  plot  the 
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time  dependence  of  intensity  for  this  pumping  level,  but  with  different  initial  condition  in  Figure 
3,  which  approaches  a  non-zero  solution. 


Figure  3:  Time-Dependence  of  Intensity  for  the  Blue  Pumping  and  Larger  Initial  Intensity 

Thus  we  can  conclude  that  the  developed  model  gives  a  correct  physical  picture  of  the 
observed  effect  and  can  serve  as  a  modeling  tool  for  designing  structures  with  desirable  bistable 
and  reconfigurable  properties.  Even  though  we  did  not  manage  to  observe  these  effects  with  the 
experimental  design  studied  during  the  Phase-1  efforts,  the  obtained  theoretical  results  proof  that 
the  general  concept  of  the  device  based  on  resonant  coupling  of  the  waveguides  is  correct  and 
can  be  realized  with  more  carefully  prepared  structure. 

2.3.  Evanescent  coupling  and  normal  modes  of  coupled  disks. 

The  second  thrust  of  the  Phase-I  efforts  was  concerned  with  observing  bistable  behavior 
in  the  system  of  coupled  microdisks.  In  order  to  provide  theoretical  supports  to  these  efforts  and 
assist  in  designing  of  the  structure  a  theoretical  model  and  relevant  software  for  analysis  of 
normal  frequencies  and  modes  of  single  and  double-disk  structures  were  developed. 

2.3.1.  Single  disks  -  resonance  frequencies  and  normal  modes 

It  is  well  known  that  microdisk  resonators  can  be  with  a  good  accuracy  modeled  as  tow¬ 
dimensional  structures,  in  which  two  different  polarizations  of  the  field  and  can  be  separated  [2], 
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As  a  result  the  problem  of  finding  resonance  frequencies  of  these  structures  is  reduced  to  solving 
a  scalar  tow-dimensional  wave  equation.  However,  when  introducing  normal  modes  of  these 
structures,  which  are  required  for  modeling  their  lasing  properties  one  needs  to  take  into  aceount 
that  miero-disk  resonators  have  relatively  low  Q-factors,  which  means  that  their  radiative  losses 
must  be  taken  seriously.  Analysis  earried  out  in  this  work  was  based  on  so  ealled  constant  flux 
modes  introduced  recently  in  Ref.  [3],  which  were  defined  in  the  following  way.  Inside  the  disks 
the  field  is  described  by  Bessel  function  Jm  (kn.r)  where  k  =  com„  /  c  determines  the  normal 

frequency  comn ,  while  outside  of  the  disk,  the  field  is  given  by  an  outgoing  Hankel  function 

Hm  ( kr ),  which  depends  on  an  external  parameter  k  =  a>/ c ,  where  spectral  parameter  does 

not  coincide  with  the  normal  frequency.  The  main  advantage  of  the  modes  defined  this  way 
compared  to  more  traditional  quasimodes  [4]  is  that,  unlike  the  latter,  they  do  not  diverge  at 
infinity  and  ean,  therefore,  be  direetly  used  to  calculate  the  outgoing  intensity  of  respective 
lasing  structures.  Normal  frequencies  for  TE  polarized  modes  in  this  case  are  defined  by 
equation: 

A,*)  I  kmn  // „  ( kR ) 

J'mAkmnndR)  nd  k  K{kR) 

where  nd  is  the  effective  refractive  index  of  the  disk.  These  frequencies  arc  eomplex  -valued 
and  depend  on  the  external  spectral  parameter.  The  real  part  of  the  normal  mode  is  defined  as 
solution  of  the  self-eonsistent  equation  Ref kmn (kr  )]  =  kr ,  while  its  imaginary  part  is  given 

as Im[&mn(&r)].  In  order  to  numerically  evaluate  these  quantities  we  developed  computer  codes 

(Appendix  A)  based  on  the  Mathematica  platform  using  its  standard  library  functions.  The 
ehoice  of  Mathematica  was  due  to  the  fact  that  it  allowed  carrying  out  calculations  of  Bessel 
functions  with  arbitrary  precision  which  was  crucial  for  finding  normal  frequencies  of  modes 
with  high  values  of  azimuthal  numberm  .  Figure  4  below  presents  the  found  frequencies  and 
their  imaginary  parts  for  modes  with  azimuthal  numbers  from  10  to  100  and  radial  numbers.  We 
also  compared  the  results  found  for  the  constant  flux  modes  with  those  found  for  quasi-modes. 
The  real  and  imaginary  parts  of  the  latter  eoineide  with  those  of  scattering  resonances  found  as 
poles  of  the  respective  scattering  amplitudes. 
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Figure  4:  Real  and  Imaginary  Parts  of  Normal  Frequencies  of  Constant  Flux  Modes  for  Azimuthal 
Mode  Numbers  from  10  to  100.  Different  Colors  Correspond  to  Different  Radial  Numbers 


The  comparison,  shown  in  Figure  5,  indicates  that  while  this  difference  might  be 
significant  under  some  circumstances,  it  is  small  enough  to  be  neglected  for  the  purpose  of 
analyzing  effects  of  evanescent  coupling  on  normal  modes  of  multiple  disks.  Since  numerical 
analysis  of  the  scattering  resonances  is  less  numerically  expansive,  the  problem  of  normal  modes 
of  the  coupled  disks  was  considered  as  the  scattering  problem. 


i 


i 


20, 
1  5- 
1  O-i 
05- 


■  1  st  radial  mode 
•  2nd  radial  mode 
A  3rd  radial  mode 


If" 

16- 


t 

14  - 

a 

12- 

<TJ 

10- 

8 

a: 

0) 

6- 

> 

4  - 

i 

% 

2  - 

a : 

0- 

-2- 

0 


X 


t 


10  20  30  40  50  60  70  80  90  100  TlO 


Angular  mode# 


Figure  5:  Relative  Differences  between  Real  and  Imaginary  Parts  of  Constant  Flux  Normal 
Frequencies  and  Scattering  Resonances:  It  is  Seen  that  for  High  Order  Modes  the  Difference 

Becomes  Very  Small 
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2.3.2.  Normal  modes  of  the  double-disk  structures 


To  model  a  double-disk  structure  it  was  found  convenient  to  solve  a  more  general 
multiple-disk  scattering  problem.  Each  disk  constituting  the  structure  was  characterized  by  its 
internal  and  scattered  fields  presented  in  the  form  of  the  following  modal  expansion 


(16) 

m 


where  /  enumerates  disks  characterized  by  radii  R  and  eenters  positioned  at  r  =  r  .  In  addition  to 

these  two  fields  one  also  needs  to  introduce  an  incident  field,  which  represents  the  initial  incident 
field  (first  term  in  the  equation  below)  and  the  fields  scattered  by  all  other  disks  (the  second 
term): 


B\nc  =  Z  a»e‘mm  JmM|r-r,|)  +  ZZ  hln  |r  "  |)  (  1  7) 


j*i  m 


Here  polar  angles  (p.  and  radial  coordinates  are  defined  in  local  coordinate  systems  associated 


with  each  7  —  //?  disk.  For  the  scattering  problem  considered  here  arguments  of  all  Bessel 
functions  eontain  the  frequency  of  the  incident  field.  In  order  to  be  able  to  use  boundary 
conditions  at  the  rim  of  the  /  —  //7  disk  the  scattered  field  of  all  disks  needs  to  be  rewritten  in  the 
coordinate  system  centered  at  the  /  —  //7  disk.  This  is  achieved  with  the  help  of  the  Graffs  formula 
for  the  Hankel  function  [5] 


(w/lr-rj)  =  Z  Hn-m  (n0kRj, )e‘[m~"KJn  (n/|r  -  r,|)< 


(18) 


where  R  and  0  are  radial  and  polar  coordinates  of  the  i  —  th  disk  in  the  coordinate  system 
centered  at  j  -th  disk.  Using  this  theorem  a  single  disk  results  is  generalized  into  the  multiple 
disk  form 


)  =  <**(<») 


+ZZ#',W/—  ("o  kRr)e‘ 


(19) 


where  a ^  (m)  is  a  single  disk  scattering  coefficients 
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In  the  ease  of  a  double-disk  system  Equation  (19)  simplifies  into 


(20) 


w 


(21) 


where  D  is  the  distance  between  the  disks’  centers,  and  it  was  assumed  that  only  the  A/-th  mode 
of  the  first  disk  is  being  excited  by  an  original  incident  field.  This  is  the  system  of  infinite 
number  of  equations,  whieh  describes  coupling  between  different  azimuthal  modes.  These 
equations  were  analyzed  under  assumption  that  all  disks  are  identical,  whieh  is  a  very  good 
approximation  for  the  structures  under  consideration.  To  solve  these  equations  it  should  be  first 
notieed  that  the  strongest  coupling  takes  plaee  between  modes  with  equal  or  elose  single  disk 
resonances.  Sinee  in  mierodisks  modes  with  azimuthal  numbers  of  the  same  magnitude  but 
opposite  signs  are  degenerate,  one  first  needs  to  consider  coupling  between  modes  described  by 
coefficients b±]M  .  Also  taking  into  aeeount  that  the  Hankel  functions  grow  fast  with  increasing 

order,  one  realizes  that  coupling  between  coefficients  b[\]  and  b{\]n  characterized  by  coupling 
parameter  H2M(kD ),  is  mueh  stronger  than  the  coupling  between  coefficients  b^  and  b(2) , 
whieh  is  characterized  by  H0(kD) .  Thus,  one  can  solve  Equation  (21)  iteratively,  taking  first  into 
aeeount  the  coupling  between  main  modes,  and  incorporating  the  rest  using  perturbation  theory 
with  coupling  coefficients  HM_n  /  H2M  as  small  parameters.  The  zero  order  approximation  takes 

into  aeeount  only  coupling  between  counter-propagating  resonant  modes  with  following 
analytical  result  for  new  resonances  modified  by  coupling: 

[aM{k)]x  =±H2M{kD)  (22) 

The  main  feature  of  this  approximation  is  that  it  preserves  degeneraey  between  two  modes,  one 
of  whieh  ean  be  described  as  symmetric  and  characterized  by  b{^+b{2^  combination  of  the 

coefficients,  while  the  other  one  is  antisymmetric  characterized  by  b{^  -b{2^  combination.  This 
equation  gives  a  simple  method  to  estimate  the  effect  of  coupling  on  normal  modes  and  was  used 
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in  conjunction  with  experimental  results  to  verify  the  effieieney  of  coupling  in  the  fabricated 
structures. 

This  approximation  was  further  improved  by  taking  into  aeeount  coupling  to  the  non¬ 
resonant  modes  in  the  first  order  of  the  perturbation  theory.  The  results  can  be  presented  in  the 
following  form 

[«„(*)]'' -±WJV(*D)+  2>„ [//„„(*£>)]=  ±  y  „(*D)  (23) 

n*  M  n±-M 

which  shows  that  now  there  arc  four  distinet  resonanee  frequencies  so  that  the  degeneracy 
between  symmetric  and  antisymmetric  modes  is  now  removed.  Using  this  result  one  ean 
determine  if  the  structure  under  consideration  requires  going  beyond  the  resonant  approximation 
for  its  modeling.  Since  in  our  experiments  only  two  frequencies  have  been  observed  we  ean 
conclude  that  interaction  with  non-resonant  modes  is  not  important  and  ean  be  neglected. 

2.4.  Summary 

Theoretical  and  modeling  efforts  with  this  Phase  I  projeet  produced  following  results,  whieh  arc 
erueial  for  the  successful  completion  of  the  next  phase  of  this  project: 

1 .  We  ascertained  the  physical  mechanism  responsible  for  multistablc  lasing  in  the  system 
of  two  waveguides  coupled  through  a  resonant  element  (ring  or  disk  resonator) 

2.  The  developed  model  of  lasing  in  this  type  of  structures  allowed  predicting  emission 
properties  of  the  structure  for  given  morphological  characteristics  and  tailor  the  structural 
parameters  to  achieve  desirable  multistable  properties 

3.  For  the  system  of  eoupled  disks  a  theoretical  model  and  necessary  eomputcr  codes  were 
developed  for  simulating  properties  of  normal  modes  of  these  structures.  The  obtained 
results  arc  crucial  for  further  modeling  of  lasing  from  these  structures. 

3.0  EXPERIMENTAL 

The  experimental  work  performed  under  this  Phase  I  grant  can  be  divided  into  two  main 
subtasks: 

1 .  The  microring  integrated  with  an  active  waveguide  system  and 

2.  The  coupled  active  microdisk  system 
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Both  these  systems  were  investigated  for  their  bi-stable  output  eharaeteristies  in  addition 
to  their  linear  transmission  and  luminescence  properties. 

3.1.  Microring  integrated  with  an  active  waveguide 

As  part  of  this  task,  we  have  been  involved  in  the  design  of  the  all  optical  flip-flop.  The 
schematic  drawing  of  the  proposed  architecture  is  shown  in  Figure  6.  Here  two  gain  media  are 
integrated  with  a  single  mieroring  resonator  as  shown  in  the  schematic  drawing.  Specific 
subtasks  investigated  are  briefly  described  below. 


Figure  6:  Schematic  Drawing  of  the  Proposed  All-optical  flip  flop  which  uses  a  passive  microring 
resonator  integrated  with  active  elements.  Also  shown  is  the  crossesction  of  the  SU8  waveguides 
and  the  glass  substrate.  The  waveguides  are  2  pmx2pm  in  si/e 

3.1.1.  Design  of  ring  resonator 

Firstly  we  investigated  the  bend  loss  mechanism  in  a  waveguide  made  using  SU8  -  a 
negative  photosensitive  polymer.  Simulations  were  performed  using  BEAMPROP,  a 
commercially  available  waveguide  modeling  tool.  Shown  in  Figure  7  is  calculated  loss  of  a  SU8 
waveguide  (2  pm  x  2  pm)  at  different  wavelengths  for  a  180  degree  bend. 
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Figure  7  Calculated  Bending  Loss  for  Different  Input  Wavelengths  as  a  Function  Of  Bend  Radius 

for  a  180  Degree  Bend 


The  blue  line  in  the  plots  above  indicates  the  power  in  the  fundamental  mode  as  a  function  of 
bending  radius  for  different  input  wavelength  values.  The  optimal  values  of  radii  chosen  for  the 
different  wavelengths  were  205,  130,  and  100  pm. 

3.1.2.  Design  of  xMultimode  Interference  Coupler. 

We  simulated  the  performance  of  the  coupler  which  will  be  used  for  coupling  the  light 
from  the  straight  waveguides  to  the  ring  resonator  using  beam  propagation  software.  Results  of 
simulations  are  shown  in  Figure  8. 


650  nm  850  nm  1550  nm 


Figure  8:  Simulations  of  Multi-mode  Interference  Couplers  with  50-50  Splitting  Ratio  for 

Different  Input  Wavelengths 


The  couplers  were  based  on  multi-mode  interference  (MMI)  scheme  and  were  designed  for  50- 
50  power  splitting  ratio. 
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3.1.3.  Resonant  frequencies  of  the  ring  resonator  with  MM1  couplers 


Following  the  design  of  the  MM1  coupler  and  the  calculations  of  optimal  bend  radius,  we 
simulated  the  resonant  frequencies  of  the  ring  resonator.  To  this  end  we  calculated  transmission 
coefficient  of  a  coupler,  T ,  using  a  model,  in  which  the  coupling  to  the  resonator  is  described  by 
a  2x2  transfer  matrix,  connecting  fields  incident  at  the  coupler  from  the  waveguide,  and  from 
the  ring.  The  coupler  is  characterized  by  power  coupling  parameters  k  ,  and  the  loss 
parameters  y ,  and  a  ,  where  former  characterizes  the  fraction  of  power  lost  at  the  coupler,  while 


the  latter  takes  into  account  absorption  in  the  resonator.  This  approach,  which  is  standard  for 
MM1  couplers  and  is  used  in  a  number  of  works,  results  in  the  following  expression  for  the 
power  transmission  coefficients  defined  as  a  fraction  of  energy  propagating  in  the  waveguide 
passing  the  coupler: 


r-O-r) 


2  —  (! -r)‘  "  exP 

'  1  ^ 

-  aL-i<X> 
v  2  J 

1  —  (1  —  at)'  ’  (l  "  exp 

1 

—  aL-i cb 

2  J 

(24) 


where  L  is  the  circumference  of  the  resonator,  and  O  =  prL ,  where  pr  is  the  propagation 

constant  in  the  ring.  Figure  9  shows  the  transmission  eharaetcristies  of  the  ring  resonator  at 
different  wavelengths.  The  minima  in  transmission  correspond  to  the  resonance  coupling  of  light 
to  the  resonator. 


650  nm  850  nm 


1550  nm 


Wavelength 


Figure  9:  Simulations  of  Resonant  Frequencies  of  The  Ring  Resonators  Designed  for  Three 
Different  Wavelength  Ranges  :  650,  850,  and  1550  nm 
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3.1.4.  Fabrication  of  microring  structure 

Microring  structures  were  realized  using  both  standard  photolithography  as  well  as 
electron  beam  lithography  following  the  design  parameters  established  above.  One  of  the  key 
components  in  this  architecture  is  the  integration  of  active  waveguide  consisting  of  colloidal 
quantum  dots  (CdSe/ZnS)  embedded  in  a  polymer  host  with  a  passive  waveguide.  We  have 
successfully  implemented  this  step  using  a  multilayer  electron  beam  lithography  proeess  with 
appropriate  alignment  marks.  Scanning  electron  microscope  image  of  an  aetive  waveguide 
integrated  with  a  passive  waveguide  is  shown  in  Figure  10  (a).  Also  shown  in  the  same  figure  is 
the  scanning  electron  microscope  image  of  the  mieroring  resonator  (b). 


Figure  10:  Scanning  Electron  Microscope  images  Of  (A)  Active  Waveguide  Integrated  With  Passive 
W  aveguide,  And  (B)  Microring  Resonator  Coupled  To  Passive  Waveguides.  Here,  The  Active 
Waveguides  Consist  Of  Colloidal  Cdse  Qds  Embedded  In  A  SU8  Matrix. 


The  aetive  passive  integrated  waveguide  structure  was  eharaetcrized  by  optical  pumping  of  the 
device  using  an  Argon  ion  laser  (488  nm).  During  the  initial  phase  of  device  fabrication,  one  of 
the  major  issues  we  had  was  that  the  colloidal  quantum  dots  were  getting  deposited  all  over  the 
wafer  despite  removing  the  polymer  using  the  developer.  This  issue  was  overcome  and  now  we 
have  been  able  to  realize  devices  where  only  the  aetive  region  has  the  quantum  dots.  Shown  in 
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Figure  11  is  an  optical  mieroseope  image  of  an  aetive  waveguide  integrated  with  a  passive 
waveguide  along  with  the  image  under  optieal  excitation  where  the  red  emission  is  observed  only 
from  the  aetive  region. 


Figure  11:  Optical  Microscope  Image  of  the  Active  Waveguide  Integrated  with  a  Passive 
Waveguide  (a)  and  Image  of  Emission  from  the  Active  W  aveguide  under  Optical  Excitation  (b) 


The  ring  resonator  was  eharaeterized  for  its  transmission  properties  using  a  broadband  source 
operating  in  the  1550  nm  wavelength  range.  Although  the  active  material  we  had  chosen  was 
CdSe  quantum  dots  that  emit  in  the  red  (620  nm),  the  ease  of  characterizing  in  the  infrared 
prompted  us  to  do  the  transmission  measurements  in  the  infrared.  Shown  in  Figure  12  is  one  of 
the  transmission  resonances  of  the  ring  resonator  in  the  1550  nm  wavelength  range. 
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Figure  12:  Transmission  Characteristics  of  the  Ring  Resonator  Shown  Indicating  Extinction  of  ~20 

dB  and  0.07  nm  FWHM 


Summary:  We  have  successfully  fabricated  a  microring  resonator  integrated  with  two  active 
waveguides.  The  integration  of  the  active  waveguide  with  the  passive  waveguide  was  shown  to 
be  sueeessful  based  on  the  emission  eharaeteristies.  Although  we  have  not  yet  observed  bistable 
lasing  behavior  from  this  integrated  device,  experiments  are  currently  underway  to  demonstrate 
this  in  the  near  future.  The  outstanding  issues  include  low  reflectivity  from  the  waveguide  facets, 
and  losses  in  the  ring  resonator. 

3.2.  Coupled  Active  Microdisks 

In  this  approach  towards  realizing  bistable  lasing,  we  use  active  microdisks  consisting  of 
CdSe/ZnS  colloidal  uantum  dots  embedded  in  a  photosensitive  polymer.  In  the  present 
demonstration,  the  polymer  used  was  SU8  (Microchem  Corp.),  a  negative  photoresist.  Patterning 
of  the  microdisks  were  achieved  using  both  soft-lithography  and  electron  beam  lithography.  Due 
to  the  resolution  limit  of  soft-lithography,  it  was  not  possible  to  realize  efficient  coupled 
microdisks  using  this  technique.  Hence  we  adopted  electron  beam  lithography  to  realize  this 
structure.  Scanning  electron  microscope  image  of  the  coupled  microdisks  embedded  with  CdSc 
quantum  dots  is  shown  in  Figure  13. 
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Figure  13:  Scanning  electron  microscope  image  of  coupled  microdisks  fabricated  using  electron 
beam  lithography.  The  disks  are  40  pm  in  diameter  and  are  separated  by  300  nm 


Here  the  separation  between  the  disks  was  varied  from  250  to  350  nm  to  tune  the  coupling 
between  them.  These  coupled  microdisks  were  characterized  for  their  emission  properties. 
Shown  in  Figure  14  is  the  image  of  emission  observed  using  a  CCD  camera  from  a  three  coupled 
microdisk  structure  along  with  the  optical  microscope  image  of  the  sample  under  investigation. 


Figure  14:  Optical  Microscope  Image  of  The  Coupled  Microdisk  Structure  Along  with  the 
Emission  from  the  Three  Coupled  Microdisk  Sample 
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We  also  performed  spectral  measurements  on  the  coupled  microdisks.  The  experimental  set  up 
used  to  characterize  these  microdisks  is  shown  in  Figure  15. 


Figure  15:  The  experimental  setup.  1  -  laser,  2  -  mirror,  3  -  microscope,  4  -  sample,  5  -  31) 
stage  with  rotation,  6  -  f=  25  mm  lens,  7  -  f=  15  mm  lens,  8  -  filter,  9  -  optical  fiber,  10  -  HR 
4000  CG-UV-NIR  spectrometer,  11  -computer. 


The  spectral  property  of  the  coupled  microdisks  under  optical  excitation  using  an  Argon- 
ion  laser  was  studied  using  a  high  resolution  spectrometer.  Optical  spectrum  observed  from  a 
two  coupled  microdisk  sample  is  shown  in  Figure  16.  Two  clear  peaks  are  observed  above  the 
spontaneous  emission  spectrum.  While  we  have  not  yet  confirmed  if  these  peaks  are  indeed 
lasing  peaks,  the  significantly  narrower  linewidths  indicate  that  there  is  gain  occurring  in  these 
systems.  In  addition,  bright  emission  observed  from  these  disks  further  substantiates  our  claims. 
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Wavelength,  nm 


Figure  16:  Emission  from  a  Two  Coupled  Microdisk  Structure  Showing  Two  Significantly  Narrow 

Peaks 


Summary:  We  have  successfully  demonstrated  coupled  microdisks  using  colloidal  quantum  dot 
composites.  These  disks  have  demonstrated  excellent  spectral  properties  and  have  shown 
spectrally  selective  emission.  Outstanding  issues  to  be  addressed  include  demonstration  of 
optical  control  of  bistable  emission,  integration  of  passive  waveguides  to  these  active  disks,  and 
short  operational  lifetime  of  these  active  microdisks. 

4.0  CONCLUSIONS  AND  RECOMMENDATIONS 

Under  the  Phase  I  grant,  the  Hybrid  Photonics  +  Queens  College  team  have  successfully 
performed  feasibility  study  on  realizing  bistable  optical  devices  using  novel  physical 
mechanisms.  Specific  tasks  accomplished  include. 

■  Developed  the  theory  of  bistable  lasing  in  a  system  consisting  of  a  single  resonator  coupled 
to  waveguides. 

■  Developed  software  that  predicts  the  lasing  modes  in  such  systems 

■  Successfully  fabricated  mieroring  resonator  side  coupled  to  active  waveguides 

■  Successfully  integrated  colloidal  quantum  dot  based  active  waveguide  with  a  passive 
waveguide. 

■  Developed  theoretical  model  and  necessary  computer  codes  for  simulating  properties  of 
normal  modes  of  coupled  micridisk  structures.  The  obtained  results  are  crucial  for  further 
modeling  of  bistable  lasing  from  these  structures. 
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■  Fabricated  eoupled  microdisk  structures  and  demonstrated  spectrally  selective  narrowband 
multiwavelength  emission. 

Recommendations  for  future  work: 

Based  on  our  preliminary  results  of  Phase  I,  the  eoupled  active  mierodisk  resonator 
system  looks  more  promising  than  the  ring  resonator  system  to  obtain  bistable  lasing.  Theoretical 
simulations  have  elearly  indicated  the  possibility  to  obtain  bistable  lasing  in  the  former  system. 
Furthermore,  experimentally  wc  have  observed  narrow  speetral  emission  from  these  eoupled 
mierodisk  structures.  Both  theory  and  experiment  favor  the  former  system  for  realizing 

In  future  work,  we  will  use  the  eoupled  microdisk  resonator  system  to  realize  the  all- 
optical  flip  FF.  The  choice  of  this  architecture  for  realizing  the  FF  is  based  on  reeent  theoretical 
and  experimental  results  we  have  obtained.  The  ring-resonator  based  structure  was  considered  as 
a  quiek  first  step  in  realizing  a  principally  new  type  of  bistability  providing  spectral  diversity,  but 
the  progress  in  this  direction  was  hindered  by  unanticipated  difficulties  sueh  as  low  light 
intensities  in  the  passive  waveguides.  Nevertheless,  as  the  main  idea  of  this  deviec  does  not 
depend  on  speeifie  [ring-resonator]  architecture,  we  have  found  recently  that  the  device  is 
realized  easier  by  microdisk  structures.  Thus  we  will  continue  pursuing  an  all-optical  FF  based 
on  the  idea  of  the  resonance  coupling  of  active  waveguides  by  a  mierodisk  resonator.  It  is 
important  to  point  out  that  the  study  of  the  waveguide  eoupled  single-disk  structure  is,  by  any 
means,  a  neeessary  intermediate  step  toward  realizing  eoupled-disk  structures.  Thus,  while  the 
main  foeus  of  our  future  work  will  be  on  eoupled-disk  based  optieal  flip-flops,  we  anticipate 
realizing  the  principally  new  type  of  single-disk  FF  en  route  toward  the  main  goal.  The  proposed 
modifications  in  our  future  efforts  refleet  realization  that  the  development  of  a  double-disk  based 
FFs  has  lesser  risks  as  well  as  guarantees  excellent  chances  of  its  successful  realization. 
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APPENDIX  A 


CALCULATIONS  OF  NORMAL  AND  SCATTERING  RESONANCES  FOR  INDIVIDUAL  AND 

COUPLED  DISKS 

Scattering  resonances : 

nO  =  0;  nd  =  0;  (*  to  make  global,  do  not  change  it  *) 

DiskTE  [mmin_,  mmax_,  mpoints_,  kmin_,  kmax_,  kpoints_,  levels2count_]  :  = 

Module  [  {out  =  Table  [{}  ,  {6}],  xl ,  al ,  j,  ReK,  ImK,  glmK,  m,  Nstep,  startTime, 

endTime,  g,  levelsO  = levels2count ,  countLevels,  R  =  SetPrecision [R ,  myPrec] } , 

Off [General: :unfl] ; 

If  [levels2count  ==  -1,  levelsO  =  10  A5;]  ; 

xl  =  SetPrecision  [Range  [kmin  ,  kmax,  (kmax  -  kmin)  /  kpoints]  ,  myPrec]  ; 
startTime  =  AbsoluteTime [ ]  ; 

If  [stepsN  t  0  ,  Nstep  =  deltaN  /  stepsN,  Nstep  =  {0 ,  0}  ;  ]  ; 

Do  [  ( *3*) 

countLevels  = 0 ; 

Print [ ’’Working  with  order=",  m,  ”  ”,  DateString []]  ; 

nO  =  SetPrecision [nlnd[ [1] ] ,  myPrec] ; 
nd  =  SetPrecision [nlnd[ [2] ] ,  myPrec] ; 
al  =  f  Zerolraag  [xl ,  m,  R]  ; 

Do [  (*2*) 

If [countLevels  >=  levelsO, 

Break [ ] ; 

]  ; 

If  [al  [  [g]  ]  *  al  [  [g  +  1]  ]  <  0, 

Do  [  (*1*) 

If  [j  ==  0, 

nO  =  SetPrecision [nlnd [ [1] ] ,  myPrec] ; 
nd  =  SetPrecision [nlnd [ [2]  ] ,  myPrec] ; 

Zerolmag  =  FindRoot  [f  Zerolmag  [x ,  m,  R]  ,  {x,  xl  [  [g]  ]  }  ,  AccuracyGoal  -> 
myGoals ,  PrecisionGoal  ->  myGoals  /  2,  WorkingPrecision  -►  myPrec  -  1]  ; 

ReK  = Zerolmag [ [1] ] [ [2] ] ; 

ImK  =  ReK  /  2  *  nO  *  nd  *  ( -  Bessel  J  [m ,  ReK  *  nd  *  R]  * 

(BesselJ  [m  -  1 ,  ReK  ★  nO  *  R]  -  Bessel  J  [m  +  1 ,  ReK  *  nO  *  R] )  +  BesselJ  [m, 

ReK  *  nO  ★  R]  *  (BesselJ [m  -  1 ,  ReK  *  nd  *  R]  -  BesselJ  [m  +  1 ,  ReK  *  nd  *  R]  )  )  ; 
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glmK  =  1  /  (ReK  *  R)  *  (2  *  ReK  *  (nO  -  nd)  *  nd  *  R  *  Bessel  J  [m  -  1 ,  ReK  *  nd  *  R]  * 
(ReK  *  nO  *  R  *  BesselY  [m  -  1 ,  ReK  *  nO  *  R]  - 

m*  BesselY  [m,  ReK*nO*R])  +  BesselJ[m,  ReK*nd*R]  * 

(ReK  *  nO  *  R  *  ( -  ReK  *  nO  *  nd  *  R  *  BesselY  [m  -  2  ,  ReK  ★  nO  ★  R  ]  -  2*  (m  ★  nO  +  nd 
2  *  m  *  nd)  *  BesselY  [m  -  1 ,  ReK  *  nO  *  R]  )  +  (4*mA2*  (nO  -  nd)  + 
ReK  A 2  *  nO  *  (nO  -  2  ★  nd)  *  nd  *  R A  2)  *  BesselY  [m,  ReK  ★  nO  ★  R]  )  )  ; 
startRe  =  ReK; 
startlm  =  ImK  /  glmK; 

r 

If  [Nstep[  [1]  ]  ==  0  |  |  deltaN[  [1]  ]  ==  0 , 
nO  =  SetPrecision[nInd[ [1] ] ,  rnyPrec] ; 

f 

nO  =  SetPrecision [nlnd [ [ 1] ]  +  I  *  j  * Nstep [ [ 1] ] ,  myPrec] ; ] ; 

(★complex  part  of  nO  *) 

If  [Nstep  [  [2]  ]  ==  0  |  |  deltaN [  [2]  ]  ==  0  , 
nd  =  SetPrecision [nlnd [ [2] ] ,  myPrec] ; 


nd  =  SetPrecision [nlnd [ [2] ]  +  I  ★ j  *  Nstep [ [2] ] ,  myPrec] ; ] ; 

(★complex  part  of  nd  ★) 
startRe  =  X[ [1] ] ; 
startlm  =  X [ [2] ] ; 

]; 

value  =  FindMinimum[ scatter  [rek,  imk,  m,  R]  ,  {rek,  SetPrecision  [startRe , 
myPrec]},  {imk,  SetPrecision  [startlm,  myPrec]},  AccuracyGoal  -►  myGoals , 
WorkingPrecision  -►  myPrec,  Method  -►  "PrincipalAxis" ]  ; 

X  =  {value [ [2] ] [[1]] [[2]] ,  value [ [2] ] [ [2] ] [[2]]}; 

Y  = value [[ 1] ] ;  (★  absolute  function  value  ★) 

,  {j,  0,  stepsN} ] ;  {*/!*) 

countLevels  +  + ; 

out [ [1] ]  =  Append [out [ [1] ]  ,  X[ [1] ]  ★  R]  ; 
out [ [2] ]  =  Append [out [ [2] ]  ,  X [ [2] ]  ★  R]  ; 
out [ [3] ]  =  Append [out [ [3] ] ,  Y]  ; 
out[ [4] ]  =  Append [out [ [4] ] ,  m]  ; 

out [ [ 5] ]  =  Append[out[[5]],  X[ [1] ] / Abs [X [ [2] ] ] ] ; 

out[[6]]  =  Append [out [[ 6] ] ,  countLevels];  (*number  of  the  level*) 

]  ; 
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,  {g,  1,  kpoints  -  1}]  ;  (*/2*) 

,  {m,  mmin,  mmax,  mpoints}]  ;  (*/3*) 

endTime  =  AbsoluteTime [ ] ; 

Print  ["Execution  time:",  (endTime  -  startTime)  ,  "  seconds"]; 

If  [out[  [1]  ]  ==  {}, 

Print ["No  solutions  found,  increase  max  K  value  or  number  of  K  points"];]; 
out  (*  return  array  *) 

]  ; 

f  Zerolmag  [x0_,  m0_,  R_]  :=  Module  [  {kO  ,  kd,  m  =  SetPrecision  [raO ,  rayPrec]  ,  x  =  xO  ,  val}, 
kO  =  x  ★  nO  ; 
kd  =  x  *  nd ; 

val  :  =  kO  *  nd  *  Bessel  J  [m  -  1 ,  kd  *  R]  *  BesselY  [m,  kO  *  R]  +  1  /  R  *  BesselJ  [ra,  kd  ★  R]  * 
(-kO  *  R  *  nd  *  BesselY  [m  -  1 ,  kO  *  R]  +  m  *  (nd  -  nO)  ★  BesselY  [m,  kO  *  R]  )  ; 

val 

]  ; 

scatter [RK_,  IK_,  m0_,  R_]  :=  Module  [  {kO  =  SetPrecision  [  (RK  +  I  *  IK)  *  nO  ,  myPrec]  , 

kd  s  SetPrecision  [  (RK  +  I  ★  IK)  ★  nd,  myPrec]  ,  val ,  m  =  SetPrecision  [mO ,  myPrec]  }  , 
val  :  =  N  [nO  *  HankelHl  [m,  kO  *  R]  *  kd  /  2  ★  (BesselJ  [m  -  1 ,  kd  *  R]  -  BesselJ  [m  +  1 ,  kd  *  R]  ) 
nd  *  BesselJ[m,  kd  ★  R]  ★  kO  /  2  * 

(HankelHl  [m  -  1 ,  kO  *  R]  -  HankelHl  [m  +  1 ,  kO  *  R]  )  ,  myPrec]  ; 

Sqrt  [Re  [val ]  A2  +  Im[val]  A2] 

]; 

(*  supplementary  functions  *) 

myFileName  :  =  "scat.  "  <>  ToString [nlnd [  [1]  ]  ]  <>"("<>  ToString [deltaN[  [1]  ]  ]  <>  " )_"  <> 
ToString  [nlnd  [  [2]  ]  ]  <>"("<>  ToString  [del  taN  [  [2]  ]  ]  <>")_"  <>  ToString  [R]  ; 
plotResults  [res_,  box_,  saveFormat_,  useColor_]  :=  Module  [  {plStyle ,  gr,  fname, 
t  =  Table [ { } ,  { 4} ] ,  styles  =  Table [ { } ,  {4} ] ,  modes ,  temp ,  k ,  color} , 

If  [res  [  [1  ]  ]  ==  {}  ,  Return [ ]  ;  ]  ; 

( *Of f [Show : : shx] ; ★)  (★  hide  warning  for  empty  graph  *) 
plStyle  =  {PlotRange  ->  All,  PlotMarkers  -►  ,  Frame  True, 

ImageSize -»  {400,  300},  Labels tyle ->  {FontSize  ->  14 } ,  Axes-»False}; 
styles [[1]]  =  {plStyle,  FrameLabel {"Re(kR)",  "Mode"}}; 
styles  [[2]]  =  {plStyle,  FrameLabel -♦  {"Re(kR)",  "Im(kR),,}}; 
styles [[3]]  =  {plStyle,  FrameLabel  -*  {"Re(kR)",  "Q-factor"}}; 
styles [[4]]  =  {plStyle,  FrameLabel  -»  {"Re(kR)",  "Absolute  Value"}}; 

If  [  !  useColor , 
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t [ [1] ]  =  ListPlot [Transpose [ {res [ [1]  ]  ,  res [ [4] ] } ] ,  styles [ [1] ] ] ; 
t  [  [2] ]  =  ListPlot [Transpose [ {res [ [1] ] ,  res [ [2] ] } ] ,  styles [ [2] ] ] ; 
t [ [3] ]  =  ListLogPlot [Transpose [ {res [ [1] ] ,  res [ [5] ]  }  ]  ,  styles [ [3] ] ] ; 
t  [  [4] ]  =  ListLogPlot [Transpose [ {res [ [1] ] ,  res [ [3] ] } ] ,  styles [ [4 ] ] ]  ; 

,  (*else*) 

modes  =  Max [res [ [ 6]  ]  ] ; 

Do  [ 

temp  =  Transpose [Select [Transpose [res] ,  U [ [6] ]  ==  k  &]  ]  ; 
color  = ColorData [ "VisibleSpectrum" ] [380 +  (750  -  380)  *  (k  -  1)  /  modes] 
t [ [1] ]  =  Show [ t [ [1]  ]  ,  ListPlot [Transpose [ {temp [ [ 1] ] ,  temp [ [4 ] ] } ] , 
{styles  [[!]],  PlotStyle  ->  color}  ]  ,  PlotRange  ->  All]  ; 

t [  [2 ] ]  =  Show [ t [ [2] ] ,  ListPlot [Transpose [ {temp [ [1] ] ,  temp [ [2] ] } ] , 
{styles  [ [2]  ]  ,  PlotStyle  -♦  color}  ]  ,  PlotRange  -♦  All]  ; 

t [  [3] ]  =  Show [ t [ [3] ] ,  ListLogPlot [Transpose [ { temp [ [1] ]  ,  temp [ [5] ] } ] 
{styles  [ [3]  ]  ,  PlotStyle  ->  color}  ]  , 

PlotRange  ->  All,  FrameTicks  ->  {Automatic,  Automatic}]  ; 

t [ [4 ] ]  =  Show [ t [ [4 ] ] ,  ListLogPlot [Transpose [ {temp [ [1] ] ,  temp [ [3] ] } ] 

{ styles  [  [4  ]  ]  ,  PlotStyle  ->  color}  ]  , 

PlotRange  -♦  All,  FrameTicks  -♦  {Automatic,  Automatic}]  ; 

,  {k,  modes} ] ; 

]  ; 

If [box  |  |  StringQ [saveFormat] , 
gr  =  GraphicsGrid [ { { t [ [1] ]  ,  t [ [2] ] } ,  { t [ [3] ] ,  t [ [4]  ]  } }  ]  ;  , 
gr  =  Graphic sColumn [ { t [ [1] ] ,  t [ [2] ] ,  t [ [3] ] ,  t [ [ 4 ] ] } ] ; 

]  ; 

If [StringQ [saveFormat] , 
fname  =  myFileName  <>  ”  .  "  <>  saveFormat; 

Print [ "Saving  graphics  to  "  <>  fname] ; 

Export  [No tebookDi rectory  [  ]  <>  fname,  gr] 

(*  "Image"  is  not  required  and  produces  strange  image  file, 

but  there  is  a  bug  in  Mathematica  7  -  it  can  not  export  ListLogPlot  *) 

]  ; 

Print [gr] ; 

]  ; 

(*  saving  data  to  the  file  *) 

saveData [da ta_,  lim_]  :  =  Module [ {k ,  1,  fname,  res,  uniqueVals,  saveData,  pos  ,  f } , 
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If  [data  [[1]]  ==  {},  Return  []  ;]  ; 

(*  check  same  value  *) 

uniqueVals  =  Union  [data  [  [1]  ]  ,  SameTest  ->  ((N[ttl,  lim]  ==N[tt2,  lim] )  &)  ]  ; 

If [Length [data [ [1] ] ]  i  Length [uniqueVals] , 

Print [ "There  are  the  same  solutions.  Try  to  perform  more  accurate  calculations 
(starting  points) . \nThese  values  for  Real  part  are:"]; 

Print [Commonest [N [data [ [1] ] ,  12] ] ] ; 

Print ["They  will  be  removed  from  saved  data"] ; 
saveData  =  Table [0,  {Length [data] } ,  {Length [uniqueVals] }] ; 

For[k  si,  ki  Length  [uniqueVals]  ,  k  +  +, 
pos  =  Position [N[data[ [1] ] , lim] ,  N [uniqueVals [ [k] ] , lim] ] ; 

For[l  =  1,  1  ^  Length  [data]  ,  1++, 
saveData [  [1]  ]  [  [k]  ]  =data[[l]]  [[pos[[l]]  [[1]]]]]; 

]  ; 

,  (*if  . .  else*) 
saveData  =  data; 

]  ; 

res  = Transpose [saveData] ; 

fname  =  myFileName  <>  "  .  txt"  ; 

f  =  OpenWrite  [NotebookDirectory  [  ]  <>  fname  ]  ; 

Do  [ 

WriteString [f ,  String Join [ 

Riffle [ToString[FortranForm[N[tt,  lim]]]  &/@res[[i]],  " \t" ] ]  <> "\n" ]  ; 

,  {i,  1,  Length [res] }] ; 

Close [f] ; 

Print ["Data  saved  to  ",  fname, 

"  in  the  format  Re(kR),  Im(kR) ,  Abs  Value,  Mode,  Qf actor"] ; 

]; 
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(★  main  part  with  the  structure  parameters  *) 
myPrec  :=50;  (*  working  precision,  should 
be  larger  than  myGoals  in  approximately  two  times  *) 
myGoals  :=  16;  (*  required  precision  for  solution  *) 
nlnd  :=  {1,  1.5};  (*  nO  and  nd  *) 

R  :=  60;  (*  radius  *) 

stepsN  :=  10;  (*  number  of  steps  for  changing  index  of  refraction  *) 

deltaN  :=  {0,  0.00001};  (*  change  in  complex  part  of  the  index  of  refraction  *) 

(*  calling  parameters:  m_min ,  m__max,  k_min ,  k_max , 

k_points ,  number  of  levels  to  find  (-1  for  10^5  levels)  *) 

results  =  DiskTE  [5 ,  15  #  10 ,  .  1 ,  2.1,  250  ,  3]  ; 

saveData  [results ,  16];  (*  data,  number  of  digits  to  save  into  file  *) 
plotResults [results ,  True,  "gif".  False]  (*  data, 

2x2  (True)  or  1x4  (False)  output,  Export  graphics  to  " format" ( "gif " , " jpg” , 
"eps". .  check  manual  for  Export  supported  formats)  file  (False  to  disable), 
use  different  color  for  each  level  (True)  ★) 
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Normal  resonances : 


nO  =  0;  nd  =  0;  (*  to  make  global,  do  not  change  it  *) 

DiskModes  [kmin_,  kmax_,  kpoints_,  kDmin_,  kDmax_ ,  kDpoints_,  levels2find_]  :  = 

Module  [{R  =  SetPrecision [R,  myPrec]  ,  ReK ,  ImK,  glmK,  startTime,  endTime,  j,  gl ,  g2 ,  k,  xO , 
xl ,  levelsO  =  levels2find,  out  =  Table [{} ,  {6}],  m  =  SetPrecision [mode ,  myPrec], 
countLevels,  Nstep,  debug  =  False,  myBreak,  pts,  Zerolmag,  al ,  rek,  imk,  same}, 

Off [General: :unfl] ; 

If  [levels2find  ==  -1,  levelsO  =  10 A 5;]  ; 

xO  =  SetPrecision  [Range  [kmin ,  kmax,  (kmax-kmin)  /  kpoints]  ,  myPrec]  ; 

xl  =  SetPrecision  [Range  [kDmin ,  kDmax,  (kDmax  -  kDmin)  /  kDpoints]  ,  myPrec]  ; 

startTime  =  AbsoluteTime  [  ]  ; 

If  [stepsN  f  0 ,  Nstep  =  deltaN  /  stepsN,  Nstep  =  {0  ,  0}  ;  ]  ; 

Do  [  (*3*) 

countLevels  =  0 ; 

k  =  xO [ [gl] ] ;  (*  value  of  the  continuous  K  *) 

Print  [  "Working  with  k:  ",  N[k,  6],  "  "  ,  gl ,  "  out  of  ",  kpoints,  "  ",  DateString  []]  ; 
nO  =  SetPrecision [nlnd [ [1] ] ,  myPrec] ; 
nd  =  SetPrecision [nlnd [ [2] ] ,  myPrec] ; 
al  =  f Zerolmag [xl ,  m,  k] ; 

Do  [  ( *2* ) 

If [countLevels  >=  levelsO, 

Break [] ; 

]  ; 

If [al [ [g2] ]  *  al [ [g2  +  1] ]  <  0 , 
myBreak  =  False; 

Do  [  (*1*) 

If  [j  ==  0, 

nO  =  SetPrecision [nlnd [ [1] ] ,  myPrec] ; 
nd  =  SetPrecision [nlnd [ [2] ] ,  myPrec] ; 

Zerolmag  =  FindRoot  [f Zerolmag [x,  m,  k]  ,  {x,  xl  [  [g2]  ]  ,  xl  [  [g2  +  1]  ]  }  ,  Method  -> 
"Brent"  ,  AccuracyGoal  -►  myGoals  /  2,  WorkingPrecision  -+  Round  [myPrec  /  1 . 2]  ]  ; 
ReK  =  Zerolmag [ [1] ]  [  [2] ] ; 

ImK  =  l/2*n0*nd*  (kA2*  Bessel  J  [m,  ReK  ★  nd  ★  R]  * 

(Bessel  J[m  -1,  k*n0*R]  -  Bessel  J  [m+1,  k*n0*R])  +  ReK  A  2  *  Bessel  J[m, 
k  *  nO  *  R]  *  (Bessel  J  [m  +  1 ,  ReK  *  nd  *  R]  -  Bessel  J  [m  -  1 ,  ReK  *  nd  *  R]  )  )  ; 
glmK  a  1/2/R*  (kA2*nO*  ndA2  *  RA2  *  ( Bessel  J[m  -  1 ,  ReK  *  nd  ★  R]  - 
Bessel  J  [m  +  1 ,  ReK  *  nd  *  R] )  *  BesselY  [m  -  1 ,  k  *  nO  *  R]  + 

(-2  *  nd  *  (ReK  *  nO  +  k  *  m  *  nd)  *  R  *  Bessel  J  [m  -  1 ,  ReK  *  nd  *  R]  +  2  /  ReK  ★ 

(-ReK*  (m-1)  *  m  *  nO  +  k  *  mA  2  *  nd  +  ReKA  3  ★  nO  *  nd  A  2  *  R  A  2)  * 

Bessel  J  [m,  ReK  *  nd  *  R]  )  *  BesselY  [m,  k  *  nO  *  R] )  ; 
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startRe  =  ReK; 
startlm  -  ImK  /  glmK, 

r 

If [Nstepf [1]  ]  ==  0  |  |  deltaN  [  [1]  ]  =  0, 
nO  =  SetPrecision [nlnd [ [1] ] ,  myPrec] ; 

r 

nO  =  SetPrecision [nlnd [ [ 1] ]  +  I  *  j  *  Nstep [ [ 1] ] ,  myPrec] ; ] ; 

(★complex  part  of  nO  ★) 

If  [Nstep  [  [2]  ]  ==  0  |  |  deltaN  [  [2]  ]  ==  0, 
nd  =  SetPrecision [nlnd [ [2] ] ,  myPrec] ; 

r 

nd  =  SetPrecision [nlnd [ [2] ]  +  I  ★  j  ★  Nstep [ [2] ] ,  myPrec] ; ] ; 

(★complex  part  of  nd  ★) 
startRe  =  X [ [1] ] ; 
startlm  =  X[ [2] ] ; 

]  ; 

If [debug, 

Print [ "Starting  from:",  startRe,  "  and  "  ,  startlm]; 

Print  ["range  ”  ,  xl  [  [g2]  ]  ,  ”  -  ” ,  xl  [  [g2  +  1]  ]  ]  ; 
pts  =  Reap [FindMinimum [scatter [rek ,  imk,  k,  m]  , 

{rek,  SetPrecision [startRe ,  myPrec]},  {imk,  SetPrecision [startlm,  myPrec]}, 
AccuracyGoal  -►  myGoals,  WorkingPrecision  ->  myPrec, 

Methods  "PrincipalAxis” ,  StepMonitor  Sow  [  {rek ,  imk}]]]  [[2,  1]  ]  ; 

Print [value] ; 

Print [pts] ; 

] ; 

value  =  FindMinimum  [scatter  [rek,  imk,  k,  m]  , 

{rek,  SetPrecision [startRe ,  myPrec]},  {imk,  SetPrecision [startlm,  myPrec]}, 
AccuracyGoal  -+  myGoals,  WorkingPrecision  -►  myPrec,  Method  -►  "PrincipalAxis " ]  ; 

X  =  {value [[2]]  [[1]]  [  [2] ]  ,  value [ [2]]  [[2]] [[2]]}; 

Y  =  value [ [1] ] ;  (★  absolute  function  value  ★) 

If [X [ [ 1] ]  <kDmin,  Print["Point  skipped  -  solution  found  below  the  min  kra"] 
myBreak  =  True ;  Break [];] ;  (★  break  of  solution  found  below  the  range  ★) 

,  { j ,  0,  stepsN} ] ;  ( ★ /l* ) 

(★  check  whether  these  values  were  found  before  ★) 

same  =  Length  [Select  [N  [Transpose  [out]  ,  10],  (8  [  [1]  ]  ==  N[X[  [1]  ]  *R,  10]  && 

8 [  [2]  ]  ==  N[X[[2]]  ★  R ,  10]  && 8 [  [4]  ]  ==  N [k  ★  R,  10])  &]  ]  ; 

If [same  ==  0,  countLevels +  + ,  Print [ "Solution  found  before"]] ; 

If  [  !  myBreak  &&  same  ==  0 , 

out [ [1] ]  =  Append [out [ [1] ]  ,  X [ [1] ]  ★  R]  ; 
out [ [2] ]  =  Append [out [ [2] ] ,  X [ [2] ]  ★  R]  ; 
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out[ [3] ]  =  Append [ out [ [3]  ]  ,  Y]  ; 

out [ [4] ]  =  Append [out [ [4]  ]  ,  k  *  R]  ; 

out [[5]]  =  Append[out[[5]],  X[ [1] ]  /  Abs [X [ [2] ] ] ] ; 

out[ [6] ]  =  Append [out [ [6] ] ,  countLevels] ;  ( *number  of  the  level*) 


]  ; 

,  {g2,  1,  kDpoints}];  (*/2*) 

,  {gl,  1,  kpoints} ] ;  (*/3*) 

endTime  =  AbsoluteTime [  ]  ; 

Print  [  "Execution  time:",  (endTime  -  star tTime)  ,  "  seconds"]; 

If  [out[  [1]  ]  {}, 

Print ["No  solutions  found,  increase  max  K  value  or  number  of  K  points"];]; 
out  (*  return  array  *) 

]; 

f  Zerolmag  [x0_,  m__,  k_]  :=  Module  [  {val ,  x  =  xO,  R  =  SetPrecision  [R,  myPrec]  }  , 

val  :  =  SetPrecision  [-  x  A  2  *  nO  *  nd  *  BesselJ [m  -  1 ,  x  *  nd  *  R]  *  BesselY  [m,  k  *  nO  *  R]  + 
Bessel  J  [m,  x*nd*R]  ★  (kA2*nO*nd*  BesselY  [m  -  1 ,  k  *  nO  *  R]  + 
m/R*  (x*nO-k*  nd)  *  BesselY  [m,  k  ★  nO  *  R] )  ,  myPrec]  ; 

val 

]  ; 

scatter  [RK_,  IK__ ,  k_ ,  m_]  :  = 

Module [ {x  =  SetPrecision [ (RK  +  I  *  IK)  ,  myPrec] ,  val ,  R  =  SetPrecision [R,  myPrec] } , 
val  :  =  SetPrecision  [nO  *  x  *  HankelHl  [m,  k*nO*R]  *l/2*x*nd* 

(BesselJ [m  -  1 ,  x  *  nd  *  R]  -  BesselJ  [m  +  1 ,  x  *  nd  ★  R] )  -  nd  *  k  *  BesselJ[m,  x  *  nd  *  R]  * 
1  /  2  *  k  *  nO  *  (HankelHl  [m  -  1 ,  k  *  nO  *  R]  -  HankelHl  [m  +  1 ,  k  *  nO  *  R]  )  ,  myPrec]  ; 

Sqrt  [Re  [val ]  A2  +  Im[val]  A2] 

]  ; 

(*  supplementary  functions  *) 

myFileName  :=  "eigen ."<>  ToString [nlnd  [  [1  ]]  ]  <> 

"  ("  <>  ToString [deltaNf [1] ] ]  <> ")_"  <>  ToString [nlnd [ [2] ] ]  <> "  ("  <> 

ToString  [deltaN[  [2]  ]  ]  <>  ")_"  <>  ToString  [R]  <>  "-M"  <>  ToString  [mode]  ; 
plotResults  [res_,  box_,  saveFormat_,  useColor_]  :=  Module  [  {pi  Style ,  gr, 

fname ,  t  =  Table [ { } ,  {4 } ] ,  styles  =  Table [ { } ,  { 4 }  ]  ,  modes ,  temp ,  k ,  color} , 

If  [res  [  [1]  ]  ==  {},  Return!];]; 

( *0f f [Show : : shx] ; *)  (*  hide  warning  for  empty  graph  *) 
plStyle  =  {PlotRange  -►  All,  PlotMarkers  -+  ,  Frame  ->  True, 

ImageSize  -►  {400,  300},  LabelStyle  ->  {FontSize  -►  14}  ,  Axes  ->  False}  ; 
styles [[1]]  =  {plStyle,  FrameLabel {"kR” ,  "Re(kraR)"}}; 
styles [[2]]  =  {plStyle,  FrameLabel {"kR",  "Irn(kmR)"}}; 
styles [[3]]  =  {plStyle,  FrameLabel  -»  {"kR",  "Q-factor"}}; 
styles [[4]]  =  {plStyle,  FrameLabel  -►  {"kR",  "Absolute  Value"}}; 

If [ l  useColor , 
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t[ [1] ]  =  ListPlot [Transpose [ {res [ [4] ] ,  res [ [1] ] }] ,  styles [ [1] ] ] ; 
t [ [ 2 ] ]  = ListPlot [Transpose [ {res [ [4] ] ,  res [[2]]}],  styles [ [2] ]] ; 
t [  [3] ]  =  ListLogPlot [Transpose [ {res [ [4] ] ,  res [ [5] ] } ] ,  styles [ [3] ] ] ; 
t [  [4] ]  =  ListLogPlot [Transpose [ {res [ [4] ] ,  res [ [3] ] } ]  ,  styles [ [4] ] ] ; 

,  (*else*) 

modes  =  Max [res [ [6] ] ] ; 

Do  [ 

temp  =  Transpose  [Select  [Transpose  [res]  ,  8  [  [6]  ]  ==  k  &]  ]  ; 

color  =  ColorData  ["VisibleSpectrum"]  [380  +  (750  -  380)  *  (k-1)  /  modes]/ 

t[ [1] ]  =  Show[ t [ [1] ] ,  ListPlot [Transpose [ {temp [ [4] ] ,  temp[ [1] ] }] , 

{styles  [  [1]  ]  ,  PlotStyle  color}]  ,  PlotRange  -►  All]  ; 

t [ [2] ]  =  Showft [ [2] ] ,  ListPlot [Transpose [ {temp [ [4] ] ,  temp [ [2] ] } ] , 

{ styles  [  [2]  ]  ,  PlotStyle  -►  color}]  ,  PlotRange  -►  All]  ; 

t [ [3] ]  =  Showft [ [3]  ]  ,  ListLogPlot [Transpose [{ temp [ [4] ] ,  tempf [5] ] }] , 
{styles  [  [3]  ]  ,  PlotStyle  -►  color}  ]  , 

PlotRange  ->  All,  FrameTicks  -►  {Automatic,  Automatic}]  ; 

t [ [4 ] ]  =  Showft [ [4] ] ,  ListLogPlot [Transpose [ {temp [ [4] ] ,  temp [ [3] ] } ] , 
{styles  [  [4]  ]  ,  PlotStyle  ->  color}]  , 

PlotRange  -►  All,  FrameTicks  ->  {Automatic,  Automatic}]  ; 

,  {k,  modes}]; 

] ; 

If [box  |  |  StringQ [saveFormat] , 
gr  =  GraphicsGrid [ { { t [ [1 ] ] ,  t [ [2] ] } ,  {t[[3]],  t[ [4] ]}}];, 
gr  =  Graphic sColumn [  { t  [  [  1  ]  ]  ,  t [ [2] ] ,  t[[3]],  t [ [ 4 ] ] } ] ; 

]  ; 

If [StringQ [saveFormat] , 
fname  =  myFileName  <>  "  .  "  <>  saveFormat; 

Print  [  "Saving  graphics  to  "<>  fname]  ; 

Export [NotebookDirectory [ ]  ofnaLme,  gr] 

(*  "Image"  is  not  required  and  produces  strange  image  file, 

but  there  is  a  bug  in  Mathematica  7  -  it  can  not  export  ListLogPlot  *) 

] ; 

Print [gr] ; 

] ; 

(*  saving  data  to  the  file  *) 

saveData  [data_,  lim_]  :=Module[{f,  k,  1,  fname,  res,  uniqueVals ,  saveData,  pos}  , 

If  [data[  [1]  ]  ==  {}  ,  Return [ ]  ;]  ; 

(*  check  saLme  value  *) 

uniqueVals  =  Union  [data  [  [1  ]]  ,  SameTest -►  ((N[ttl,  lim]  ==N[82,  lim] )  &)]; 

If [Length [data [ [1] ] ]  ^  Length [uniqueVals] , 

Print ["There  are  dublicated  solutions.  Try  to  perform  more  accurate  calculations 
(starting  points) . \nThese  values  for  Real  part  are:"]; 
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Print [Commonest [N [data [ [1] ] , 12] ] ] ; 

Print["They  will  be  removed  from  saved  data"]; 
saveData  =  Table [0,  {Length [data] } ,  {Length [unique Vais] }] ; 

For[k  =  1,  k  £  Length [ unique Vais] ,  k  +  +, 
pos  =  Position [N [data [ [1] ] ,  lim]  ,  N [uniqueVals [ [k] ] ,  lim] ] ; 

For[l  =  1,  1  £  Length  [data]  ,1++, 
saveData  [  [1]  ]  [  [k]  ]  =data[[l]]  [  [pos  [  [1]  ]  [  [1]  ]  ]  ]  ]  ; 

]  ; 

,  (*if  . . .  else*) 
saveData  =  data; 

]  ; 

(*res=Transpose [saveData] ;*) 

res  =  SortBy  [Transpose  [saveData]  ,  #[6]  &]  ;  (*  sort  by  number 

of  the  level  to  manually  select  required  levels  from  saved  data  *) 
fname  =  myFileName  <>  "  .  txt"  ; 
f  =  OpenWrite  [NotebookDirectory  [  ]  <>  fname  ]  ; 

Do  [ 

WriteString [ f , 

StringJoin [Riff le [ToString[FortranForm[N [U,  lim]]]  &/@res[[i]],  "\t" ]]  <>"\n" ] 
,  {i,  1,  Length [res] }] ; 

Close [f] ; 

Print ["Data  saved  to  ",  fname, 

"  in  the  format  Re(k_m  R) ,  Im(k_m  R) ,  Abs  Value,  k,  Qfactor" ]; 

]  ; 


(*  main  part  with  the  structure  parameters  *) 
myPrec  :=50;  (*  working  precision,  should 
be  larger  than  myGoals  in  approximately  two  times  *) 
myGoals  :=16;  (*  required  precision  for  solution  *) 
nlnd  :=  {1,  1.5};  (*  nO  and  nd  *) 

R  :  =  60 ;  (*  radius  *) 

stepsN  :  =  15 ;  (*  number  of  steps  for  changing  index  of  refraction  *) 

deltaN  :=  {0,  0.00001};  (*  change  in  complex  part  of  the  index  of  refraction  *) 

mode  :=130;  (*  bessel ' s  order  *) 

(*  calling  parameters:  k_min ,  k_max,  k_jpoints ,  km__min ,  km_max ,  km_points, 

number  of  low  levels  starting  from  km_min  (set  to  -1  to  show  all  levels  within  range)  *) 
results  =  DiskModes  [1.1,  2.1,  15,  .7,  2.5,  100,  6]  ; 

saveData  [results ,  16];  (★  data,  number  of  digits  to  save  into  file  *) 
plotResults [results ,  True,  "gif",  False]  (★  data, 

2x2  (True)  or  1x4  (False)  output,  Export  graphics  to  "format" 

("gif " , " jpg" , "eps" . .  check  manual  for  Export  supported  formats)  file  (False  to  disable), 
use  different  color  for  each  level  (True)  *) 

Do  [ 

mode  :  =  m  ; 

Print  [  "Working  with  mode  ",  N[m]  ,  "  ",  DateString  []]  ; 
results  =  DiskModes  [  .  01 ,  1.5,  50  ,  .  01 ,  1.5,  50 ,  3]; 
saveData [results ,  16]; 

plotResults [results ,  True,  False,  False]  ; 

,  {m,  55,  95,  10}] 

Working  with  mode  5.  Sun  30  Aug  2009  20:14:03 
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Coupled  disks : 

nO  =  0;  nd  =  0;  (*  to  make  global,  do  not  change  it  *) 

DiskTE  [mVal_,  kmin_ ,  kmax_,  kpoints_,  levels2find_]  :  = 

Module  {out  =  Table  [{},  {7}],  k,  d,  NN,  same,  levelsO  =  levels2find, 

klnt,  xl,  al,  j,  ReK,  ImK,  glmK,  m,  Nstep,  startTime,  endTime,  g, 
countLevels ,  R  =  SetPrecision [R ,  myPrec]  ,  oldStart  =  {{},  {}}}, 

Off [General: :unfl] ; 

If  [level s2find  ==  -  1 ,  levelsO  =  10  A  5;]  ; 

xl  =  SetPrecision  [Range  [kmin ,  kmax,  (kmax  -  kmin)  /kpoints]  ,  myPrec] 
startTime  =  AbsoluteTime [ ] ; 

If  [stepsN  i  0  ,  Nstep  =  deltaN  /  stepsN,  Nstep  =  { 0 ,  0} ;  ]  ; 
m  =  SetPrecision [mVal ,  myPrec] ; 
d  =  SetPrecision [dD ,  myPrec] ; 

NN  =  SetPrecision [Nval ,  myPrec] ; 

Do  [  (*3*) 

countLevels  =  0 ; 

Print [ "Working  with  k=",  klnt,  "  ",  Datestring []] ; 

nO  =  SetPrecision [nlnd [ [1]  ]  ,  myPrec] ; 
nd  =  SetPrecision [nlnd [ [2] ] ,  myPrec] ; 
k  =  SetPrecision [klnt ,  myPrec] ; 
al  =  f Zerolmag  [xl ,  m,  k,  R]  ; 

Do[  (*2*) 

If [countLevels  >=  levelsO, 

Break [ ] ; 

]  ; 

if [al [ [g] ]  *  al [ [g  + 1] ]  <  0, 

Do [  (*1*) 

If[j  ==  0, 

nO  =  SetPrecision [nlnd [ [1] ] ,  myPrec] ; 
nd  =  SetPrecision [nlnd [ [2] ] ,  myPrec] ; 

Zerolmag  =  FindRoot  [f Zerolmag  [x ,  m,  k,  R]  , 

{x ,  (xl  [  [g  +  1]  ]  +  xl  [  [g]  ] )  /  2}  ,  AccuracyGoal  myGoals  , 
PrecisionGoal  -►  myGoals,  WorkingPrecision  ->  myPrec  -  1]  ; 
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ReK  =  Zerolmag [ [1] ]  [  [2] ]  ; 


r  k  7T  -I 

ImK  =  -  1  +  signEq  ★  2  ★  BesselJ  [2  ★  m,  d  *  ReK  *  nO]  *  Cos  I -  ; 

1  1  +  NN  J 

glmK  =  D  ^  (  ( -  BesselJ  [  - 1  +  m,  kO  *  nd  *  R]  +  BesselJ  [1  +  m,  kO  *  nd  *  R]  )  *  BesselY  [m, 

kO  *  nO  ★  R]  +  BesselJ[m,  kO  *  nd  *  R]  *  (BesselY  [-  1  +  m,  kO  *  nO  *  R]  - 
BesselY[l  +  m,  kO  *  nO  *  R] ) )  /  (BesselJ [m,  kO  *  nd  *  R]  *  (-BesselJ[ 

- 1  +  m,  kO  *  nO  *  R]  +  BesselJ [1  +  m,  kO  *  nO  *  R] )  +  BesselJ [m,  kO  * 
nO  *  R]  *  (Bessel J [-  1  +  m,  kO  *  nd  *  R]  -  BesselJ  [1  +  m,  kO  *  nd  *  R]  )  )  + 

signEq  *  2  *  BesselY  [2  *  m,  d  *  kO]  ★  Cos  i  - 1  ,  kO  /  .  kO  ->  ReK; 

1  1  +NN  J  J 

startRe  -  ReK; 
startlm  =  -  ImK  /  glmK; 

t 

If [Nstep[ [1] ]  =  0  |  |  deltaN [ [1] ]  =  0, 
nO  =  SetPrecision[nInd[ [1] ]  ,  myPrec] ; 

t 

nO  =  SetPrecision [nlnd [ [1] ]  +  X  *  j  * Nstep [ [1] ] ,  myPrec] ; ] ; 

(*  complex  part  of  nO  *) 

If  [Nstep  [  [2]  ]  ==  0  |  |  deltaN [  [2]  ]  ==  0  , 
nd  * SetPrecision [nlnd [ [2] ] ,  myPrec] ; 

/ 

nd  =  SetPrecision [nlnd[ [2] ]  +  I  *  j  * Nstep [ [2] ] ,  myPrec] ; ] ; 

(*  complex  part  of  nd  *) 
startRe  =  X[ [1] ] ; 
startlm  =  X[ [2] ] ; 

]■• 

(★Print [ {startRe , startlm} ] ; ★) 

value  =  FindMinimum [ scatter [rek ,  imk,  m,  k,  R] ,  {rek,  SetPrecision [startRe , 
myPrec]  }  ,  {imk,  SetPrecision  [startlm,  myPrec]  }  ,  AccuracyGoal  -►  myGoals, 
WorkingPrecision  -►  myPrec,  Method  ’’PrincipalAxis"]  ; 

( *Print [value] ;★) 

X  =  {value [ [2] ] [ [1] ] [[2]] ,  value[[2]][[2]][[2]]J; 

Y  =  value[[l]]  (★  absolute  function  value  ★) 
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,  {j,  o,  s 


tepsN) |  ;  (*/l*) 


(*  check  whether  these  values  were  found  before  *) 
same  =  Length [Select [N [Transpose [out] ,  10]  , 

(«[[!]]  ==  N[X[  [1]  ]  *R,  10]  &&»[[2]]  =  N[X[[2]]  *R,  10]  &&«[[7]]  ==  klnt)  &]  ] 
lf[same==0,  countLevels ++ ,  Print  ["Solution  found  before"]]; 

out[ [1]  ]  =  Append [ out [ [1]  ]  ,  X[ [1] ]  *  R]  ; 
out[  [2]  ]  =  Append  [out  [  [2]  ]  ,  X[  [2]  ]  *R]  ; 
out [ [3] ]  =  Append [out [ [3] ] ,  Y]  ; 
out[ [4] ]  = Append[out[ [4] ] ,  m]  ; 

out [ [5]  ]  =  Append [ out [ [5]]  ,  X[ [1] ]  /  Abs[X[ [2] ]]]  ; 

out[[6]]  = Append[out[ [6] ] ,  countLevels];  (*number  of  the  level*) 
out[ [7] ]  =  Append [out [ [7] ] ,  klnt] ; 

,  {g,  1,  kpoints  -  1}]  (*/2*) 

,  [klnt,  1,  Nvaljj;  (*/3*) 
endTime  =  AbsoluteTime [ ] ; 

Print  [  "Execution  time :  "  ,  (endTime  -  star  tTime)  ,  "  seconds"]; 

If  [out[  [1]  ]  ==  {}, 

Print ["No  solutions  found,  increase  max  K  value  or  number  of  K  points"];]; 
out  (★  return  array  *) 


f Zerolmag  [x_,  m_ ,  k_ ,  R_]  :  = 

Module ^{k0,  kd,  val ,  NN  =  SetPrecision [Nval ,  myPrec] ,  d  =  SetPrecision [dD ,  myPrec] } , 

kO  =  x  *  nO  ; 
kd  =  x  *  nd ; 

val  =  (  (BesselJ  [  - 1  +  m,  kd  *  R]  -  Bessel  J  [1  +  m,  kd  *  R]  )  *  BesselY  [m,  kO  *  R]  + 

BesselJ[m,  kd  ★  R]  *  ( -BesselY  [-  1  +  m,  kO  *  R]  +  BesselY  [1  +  m,  kO  *  R]  )  )  / 

(BesselJ[m,  kd  *  R]  *  (BesselJ  [- 1  +  m,  kO  *  R]  -  BesselJ  [1  +  m,  kO  ★  R]  )  + 

BesselJ[m,  kO  *  R]  *  (-BesselJ  [-1  +  m,  kd*R]  +  BesselJ  [1  +  m,  kd*R]))  + 


signEq  *  2  *  BesselY  [2  *  m,  d  *  kO]  *  Cos 


1  'T' 


val 


A -  14 


scatter  [RK_,  IK_,  mO_,  k_,  R_]  :=  Module  [{kO  =  SetPrecision  [  (RK  +  I  ★  IK)  *  nO  ,  myPrec]  , 
kd  =  SetPrecision  [  (RK  +  I  *  IK)  *  nd ,  myPrec]  ,  val ,  m  =  SetPrecision  [mO  ,  myPrec]  , 
d  =  SetPrecision [dD ,  myPrec] ,  NN  =  SetPrecision [Nval ,  myPrec] } , 
val  =  N  [  (nO  *  HankelHl  [m,  kO  *  R]  ★  kd  /  2  *  (Bessel  J  [m  -  1 ,  kd  *  R]  -  Bessel  J  [m  +  1 ,  kd  *  R]  ) 
nd  *  BesselJ[m,  kd  ★  R]  *  kO  /  2  * 

(HankelHl  [m  -  1 ,  kO  *  R]  -  HankelHl  [m  +  1 ,  kO  *  R]  ) )  / 

(nd  ★  BesselJ [m,  kd  *  R]  *  kO  /  2  *  (Bessel J [m  -  1 ,  kO  *  R]  -  Bessel  J[m  +  1 ,  kO  *  R] )  - 
nO  *  BesselJ  [m,  kO  ★  R]  *  kd  /  2  *  (BesselJ  [m  -  1 ,  kd  *  R]  -  BesselJ  [m  +  1 ,  kd  *  R]  )  ) 
signEq  *  2  *  HankelHl  [2  ★  m,  kO  *  d]  *  Cos  [n  *  k  /  (1  +  NN)  ]  ,  myPrec]  ; 

Sqrt  [Re  [val]  A2  +  Im[val]  A2] 

]  ; 

(*  supplementary  functions  *) 

myFileName  :=  "scat .  general .  sign”  <>  ToString  [signEq]  <> 

"  .  "  <>  ToString  [nlnd  [  [1]  ]  ]  <>”("<>  ToString  [del taN  [  [1]  ]  ]  <>”)_"<> 

ToString  [nlnd  [  [2]  ]  ]  <>”("<>  ToString  [del taN  [  [2]  ]  ]  <>”)_”<>  ToString  [R]  ; 
plotResults  [res_,  box_,  saveFormat_,  useColor__]  :=  Module  [  (plStyle ,  gr,  fname, 
t  =  Table[ { } t  {4}]  ,  styles  =  Table [{} ,  {4}] ,  modes,  temp,  k,  color) , 

If  [res  [  [1]  ]  ==  {}  ,  Return []  ;]  ; 

(★Off [Show: : shx] ; ★  )  (★  hide  warning  for  empty  graph  ★) 
plStyle  =  (PlotRange  All,  PlotMarkers  ->  ,  Frame  -►  True, 

ImageSize  ->  (400 ,  300}  ,  LabelStyle  ->  (FontSize  14 )  ,  Axes  ->  False)  ; 
styles  [[1]]  =  {plStyle,  FrameLabel  -►  {"k",  "Re(kR)"}}; 
styles [[2]]  =  {plStyle,  FrameLabel ->  {"k",  "Im(kR)"}}; 
styles [[3]]  =  {plStyle,  FrameLabel ->  {”k”,  ”Q- factor ”}} ; 
styles  [[4]]  =  {plStyle,  FrameLabel  ->  {”k",  ” Absolute  Value”}}; 

If [  I  useColor , 

t [  [1] ]  =  Lis tPlot [Transpose [ {res [ [7] ] ,  res [ [1] ] } ]  ,  styles [ [1] ] ] ; 
t [ [2] ]  =  ListPlot [Transpose [ {res [ [7] ]  ,  res[[2]]}] ,  styles[[2]]]  ; 
t  [  [3] ]  =  ListLogPlot [Transpose [ {res [ [7] ] ,  res[[5]]}] ,  styles[[3]]] ; 
t  [  [4] ]  =  ListLogPlot [Transpose [ {res [ [7] ] ,  res [ [3] ] } ] ,  styles [ [4] ] ] ; 

,  (*else*) 

modes  =  Max [res [ [6]  ]  ]  ; 

Do  [ 

temp  =  Transpose  [Select  [Transpose  [res]  ,  »[[6]]  ==k&]]; 

color  =  ColorData  [  ”  VisibleSpectrum”  ]  [380  +  (750  -  380)  *  (k  -  1)  /  modes]; 

t[  [1] ]  =  Show[t[ [1] ] ,  ListPlot [Transpose [ {temp [ [7] ]  ,  temp[ [1] ] }] , 
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{ styles  [  [1]  ]  ,  Plots tyle  ->  color}]  ,  PlotRange  -►  All]  ; 

t [ [2] ]  =  Show[t [ [2] ] ,  Lis tPlot [Transpose [ {temp[ [7] ]  ,  temp [ [2] ] } ] , 
{styles  [  [2]  ]  ,  PlotStyle  ->  color}]  ,  PlotRange  -►  All]  ; 

t [ [3] ]  =  Show[t [ [3] ] ,  Lis tLogPlot [Transpose [ {temp[ [7] ]  ,  temp [ [5] ] } ]  , 

{  styles  [  [3]  ]  ,  PlotStyle  ->  color}  ]  , 

PlotRange  -►  All,  FrameTicks  -»  {Automatic,  Automatic}]  ; 

t [  [4] ]  =  Show[t [ [4] ] ,  Lis tLogPlot [Transpose [ { temp[ [7] ] ,  temp[ [3] ] } ] , 
{styles  [  [4 ]  ]  ,  PlotStyle  -+  color}  ]  , 

PlotRange  ->  All,  FrameTicks  ->  {Automatic,  Automatic}]  ; 

,  {k,  modes}]; 

]  ; 

If [box  |  |  StringQ [saveFormat] , 
gr  =  GraphicsGrid [ { { t [ [1] ] ,  t[[2]]>,  {t[ [3] ] ,  t[ [4] ]>}];, 
gr  =  Graphic sColumn [ { t [ [1] ] ,  t [ [2] ]  ,  t [ [3] ]  ,  t[[4]]}]; 

]  ; 

If [StringQ [saveFormat] , 
fname  =  myFileName  <>"."<>  saveFormat; 

Print [ "Saving  graphics  to  "<>  fname]  ; 

Export  [NotebookDirectory  [  ]  <>  fname,  gr] 

(*  "Image”  is  not  required  and  produces  strange  image  file, 

but  there  is  a  bug  in  Mathematica  7  -  it  can  not  export  Lis tLogPlot  *) 

] ; 

Print [gr] ; 

] ; 

(*  saving  data  to  the  file  ★) 

saveData  [data_,  lim_]  :  =  Module  [{k,  1,  fname,  res,  uniqueVals ,  saveData ,  pos ,  f }  , 

If  [data  [  [1]  ]  ==  {}  ,  Return  []  ;]  ; 

(*  check  same  value  *) 

uniqueVals  =  Union  [data  [  [1]  ]  ,  SameTest  -►  ((N[ttl,  lim]  ==N[82,  lim] )  &)  ]  ; 

I f [ Length [ da ta [ [ 1 ] ] ]  t  Length [ uniqueVal s ] , 

Print ["There  are  the  same  solutions.  Try  to  perform  more  accurate  calculations 
(starting  points) . \nThese  values  for  Real  part  are:"]; 

Print [Commonest [N [data [ [1] ] ,  12] ] ] ; 

Print ["They  will  be  removed  from  saved  data"]; 
saveData  =  Table [0,  {Length [data] } ,  {Length [uniqueVals] }] ; 

For[k  =  1,  k  £  Length [uniqueVals]  ,  k  +  +. 
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pos  =  Position [N [data [ [1] ] ,  lim] ,  N [uniqueVals [ [k] ] ,  lim] ] ; 

For[l  =  1,  1  £  Length  [data]  ,1++, 
saveData [ [1] ] [ [k] ]  =  data [ [1] ] [ [pos [ [1] ] [ [1] ]]]] ; 

]; 

,  (*if  . . .  else*) 
saveData  =  data; 

]  ; 

res  =  Transpose [saveData] ; 
fname  =  myFileName  <  >  "  .  txt"  ; 

f  =  OpenWrite  [NotebookDirectory  [  ]  <>  fname  ]  ; 

Do  [ 

WriteString [f ,  StringJoin [ 

Riffle  [ToString  [FortranForm[N[tt ,  lim]]]  &/@res[[i]],  "\t"  ]  ]  <>  "\n"  ]  ; 
,  {i,  1,  Length [res] }] ; 

Close [f ] ; 

Print ["Data  saved  to  ",  fname, 

"  in  the  format  Re(kR),  Im(kR) ,  Abs  Value,  Mode,  Qfactor,  k"  ]  ; 

]  ; 


(*  main  part  with  the  structure  parameters  *) 
myPrec  :  =  50 ;  (*  working  precision,  should 
be  larger  than  myGoals  in  approximately  two  times  *) 
myGoals  :=  16;  (*  required  precision  for  solution  *) 
nlnd  :=  {1,  1.5};  (*  nO  and  nd  *) 

R  :=  60;  (*  radius  *) 

dD  :=  (2  +  0.00001)  *R;  (*  distance  between  disks  *) 

Nval  :=2;  (*number  of  disks*) 

stepsN  :=  10;  (*  number  of  steps  for  changing  index  of  refraction  *) 

deltaN  :={0,  0.00001};  (*  change  in  complex  part  of  the  index  of  refraction  *) 

signEq  :=  1;  (*  defines  sign  in  front  of  t  in  equation  *) 


(★  calling  parameters:  m,  k_min,  k_max,  k_points  , 
number  of  levels  to  find  (-1  for  10*5  levels)  *) 
results  =  DiskTE  [100 ,  0.5,  4,  200,  2]; 

saveData  [results ,  16];  (*  data,  number  of  digits  to  save  into  file  *) 
plotResults [results ,  True,  "gif".  False]  (★  data, 

2x2  (True)  or  1x4  (False)  output.  Export  graphics  to  "format" 

( "gif " , " jpg" , "eps" . .  check  manual  for  Export  supported  formats)  file  (False  to  disable), 
use  different  color  for  each  level  (True)  *) 
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AppendColumn [vector_List ,  matrix_List]  :=  MapThread[ Append,  (matrix,  vector}] 
fData  =  {}  ;  distances  =  {}  ; 

Do  [ 

dD  :  =  dist; 

results  = DiskTE [100 ,  0.1,  4,  50,  1]; 
fData  =  Append [fData,  results] ; 
distances  =  Append [distances ,  dist] ; 

,  [dist,  2 . 000001  *  R,  2.1*R,  0.01*R}]; 
fData  =  AppendColumn [distances ,  fData] ; 
saveData [fData ,  16]; 


